Chronic inflammatory arthritis drives systemic changes in circadian energy metabolism

Significance Rheumatoid arthritis (RA) is a debilitating chronic inflammatory disease in which symptoms exhibit a strong time-of-day rhythmicity. RA is commonly associated with metabolic disturbance and increased incidence of diabetes and cardiovascular disease, yet the mechanisms underlying this metabolic dysregulation remain unclear. Here, we demonstrate that rhythmic inflammation drives reorganization of metabolic programs in distal liver and muscle tissues. Chronic inflammation leads to mitochondrial dysfunction and dysregulation of fatty acid metabolism, including accumulation of inflammation-associated ceramide species in a time-of-day–dependent manner. These findings reveal multiple points for therapeutic intervention centered on the circadian clock, metabolic dysregulation, and inflammatory signaling.

joint) (1). Tissue collection occurred at ZT0, 4,8,12,16 or 20 on day 7 after the development of disease symptoms in CIA mice, or at matched time points in naive control mice. Terminal blood was collected in MiniCollect K3EDTA tubes (Greiner Bio-One) and stored on ice prior to centrifugation at 3000 x g for 10 minutes at 4ºC. Plasma was transferred to a cryovial and flash frozen in liquid nitrogen.
Tissue was transferred to cryovials and flash frozen in liquid nitrogen immediately after dissection.

Lipopolysaccharide (LPS) challenge
Female C57BL6 Kmt2c flox mice (aged 12-22 weeks) were housed in light controlled cages and exposed to 12h:12h L:D cycles for three weeks prior to the experiment.

In vivo phenotyping
Body composition was measured using an EchoMRI Body Composition Analyzer E26-258-MT machine (Echo Medical Systems). Activity and body temperature measurements were made following surgical implantation of radio telemetry devices (TA-F10, Data Sciences International) i.p. 7 to 10 days after the first CIA immunisation. Mice recovered for at least 4 days before being singly housed for telemetry recording. Measurements were recorded from day 14 after initial CIA immunisation for an asymptomatic period of up to one week, until booster immunisation on day 21. Recording was then resumed upon the development of symptoms. Food intake was monitored from day 14 by weighing remaining food pellets of singly housed mice at the start and end of each light phase.
Blood samples for measurement of fasting insulin and glucose measurement were collected eight days after the development of CIA symptoms. Food was withdrawn from the mice at ZT0. Blood was collected following removal of the tail tip at ZT8.
Glucose concentration was measured immediately using an Aviva Accu-Chek meter (Roche). Remaining blood was centrifuged at 3000 x g for 10 minutes at 4ºC then stored at -80ºC for later analysis. Insulin level was measured by ELISA (Merck Millipore, EZRMI-13K) according to the manufacturer's instructions.

Histology
Paws for histological analysis were skinned, fixed overnight in formalin and then decalcified by incubation in Osteosoft Mild Decalcifier solution (VWR) for two weeks. Tissue was then paraffin embedded, sectioned and stained with either H&E (for cellular structures) or Safranin O (for cartilage) according to standard protocols. Slides were imaged by the University of Manchester Bioimaging Core Facility using a Pannoramic 250 Flash slide scanner (3DHistech) using a 20x/0.80 Plan Apochromat objective (Zeiss).

Cytokine analysis
For corticosterone measurements, serial tail blood samples were obtained over 48 hours. Blood was immediately centrifuged at 3000 x g for 10 minutes at 4ºC, then plasma was diluted 100-fold in PBS and frozen for later analysis using the Corticosterone ELISA kit (Abcam, ab108821). Terminal plasma samples were analysed using the BioPlex Pro Mouse Chemokine 33-plex panel (BioRad, catalogue reference 12002231). Plasma was diluted 1 in 4 in standard diluent buffer before mixing with assay beads. Samples were analysed on a BioPlex 200 machine. IL6 level was measured using the mouse IL6 ELISA kit (Abcam, ab100712). Plasma was diluted 1 in 4 in dilution buffer prior to application to the ELISA plate. Absorbance was measured using a GloMax Multi Detection System plate reader (Promega).

RNA extraction
Joint tissue was ground with liquid nitrogen using a pestle and mortar, then transferred to a Lysing Matrix D tube containing Trizol. Tissue was homogenised using a BeadMill homogeniser (3 x 4 m/s for 40s). RNA was extracted using chloroform then precipitated using isopropanol. After washing with 75% ethanol, the RNA pellet was resuspended in RNase-free water. RNA was purified from liver samples using the SV Total RNA kit (Promega), according to the manufacturer's instructions. Tissue was homogenised using Lysing Matrix D tubes loaded into a BeadMill homogeniser (4 m/s for 20s). RNA was eluted in RNase-free water. RNA was purified from muscle samples using the ReliaPrep Tissue kit (Promega), following the manufacturer's instructions for the purification of RNA from fibrous tissue. Muscle tissue was homogenised using the same protocol as for joint, and eluted in RNase-free water.

RNA-seq
Sequencing library preparation and sequencing was performed by the University of Manchester Genomic Technologies Core Facility. Sample quality was determined using a 2200 TapeStation (Agilent Technologies). Libraries were generated using the TruSeq Stranded mRNA assay (Illumina, Inc.) according to the manufacturer's protocol. The multiplexed libraries were analysed by pairedend sequencing on a HiSeq 4000 instrument (76 + 76 cycles, plus indices), then demultiplexed and converted using bcl2fastq software (v2.17.1.14, Illumina).

Data analysis
Differential expression analysis was run in R (2) using edgeR (v3.30.3). Genes were considered to the differentially expressed (DE) if the false discovery rate (FDR) was less than 1 x 10 -20 (for joint) or less than 0.001 (for muscle and liver). Differential rhythmicity analysis was performed using compareRhythms R package (v0.99.0, (3)). A model selection approach was used with genes being assigned to either arrhythmic, gain of rhythm, loss of rhythm, same rhythm in both, or a change in rhythm. A probability of being in a category of at least 0.6 was required for assignment. To avoid losing genes that were clearly differential rhythmic an extra category was used where the probability of either being a gain, loss or change in rhythm was greater than 0.6. Additional rhythmicity analysis was run using the JTK-cycle functionality of MetaCycle (v1.2.0), with period length fixed to 24 hours.
Genes were considered to oscillate in naïve and/or CIA if the JTK-cycle adjP < 0.05 for one or both conditions. For comparison of JTK-cycle and compareRhythms analysis of genes classified as losing or gaining rhythmicity with CIA, raw counts were normalised by subtracting the mean of each treatment (naïve or CIA), and dividing by the standard deviation across both treatments. Acrophase was calculated for each gene by fitting a sine wave (period constrained to 24 hours) to the normalised counts from the rhythmic group, and genes were aligned to this acrophase.
Pathway analysis of the gene lists defined above used the Enrichr web tool (4,5) to detect significantly enriched pathways within the WikiPathways Mouse database  Supplementary Table S2.

Phosphoenrichment
Phosphoenrichment was done on an Agilent Bravo AssayMAP robot using Fe(III)-NTA cartridges (7) with slight adaptations. Cartridges were primed in ACN with 0.1% TFA and equilibrated with 80% ACN in 0.1% TFA. Peptides were loaded onto the cartridges followed by a wash with 80% ACN 0.1% TFA. Phosphopeptides were eluted with 1% NH3 and dried down in a vacuum centrifuge.

Mass spectrometry
Peptides were injected into a liquid chromatography-mass spectrometry (LC-MS) system comprised of a Dionex Ultimate 3000 nano LC and a Thermo Fusion Lumos. Peptides were separated on a 50-cm-long EasySpray column (ES803; Thermo Fisher) with a 75-µm inner diameter and a 60 minute gradient of 2 to 35% acetonitrile in 0.1% formic acid and 5% DMSO at a flow rate of 250 nL/min. Data was acquired with the APD peak picking algorithm at a resolution of 120,000 and AGC target of 4e5 ions for a maximum injection time of 50 ms for MS1 spectra.
The most abundant peaks were fragmented after isolation with a mass window of 1.6 Th with normalized collision energy 28% (HCD). MS2 spectra were acquired in the ion trap in rapid scan mode for a maximum injection time of 35 ms.

Phosphoproteome data analysis
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD032723 (8). RAW files were processed in Maxquant. Identified phosphosites (phospho(STY).txt), were initially viewed and filtered in using the Perseus Framework. Potential contaminants and reverse peptides were removed.
Phosphosites were filtered using a localisation probability of x>0.75, log2 transformed and further filtered to remove missing values, where sites with fewer than 15 valid (not N/A) values in a group (either CIA or naïve) were excluded.
Missing values were imputed using random numbers drawn from a normal distribution with a width of 0.3 and down shift of 1.8. Ion intensities of identified phosphopeptides were normalized between each sample using trimmed means of M-values function from the edgeR (v3.30.3) R package. Differential phosphorylation analyses between groups were conducted using edgeR using a 5% FDR (9). Protein kinases analysis was performed using kinswingR package (10), using the curated kinase substrate sequences mouse dataset from PhosphoSitePlus (11). Phosphopeptide differential rhythmicity analysis was performed using compareRhythms (v0.99.0, (3)) in the same way as previously described for the RNA-seq data.

Western blotting
Mouse livers were lysed with western blot lysis buffer (NaCl 150mM (Sigma

Metabolomic analysis
Global metabolite profiling of liver, muscle and plasma samples was performed by Metabolon (Durham, NC, USA). Samples were analysed using the HD4 platform, which uses ultrahigh performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) methods for metabolite detection. Peaks were quantified by area under the curve, and normalised to set metabolite median value to one. Analysis by two-way ANOVA was used to identify metabolites showing altered level with treatment and/or time. Metabolites were considered to be significantly altered at a time point if the one-way ANOVA contrast false discovery rate (q-value) was < 0.05. Metabolomic summary data and details of statistical tests are provided in Supplementary Dataset S2. Metabolite differential rhythmicity analysis was performed using CompareRhythms (v0.99.0,(3)) in the same way as previously described for the RNA-seq data.
For inter and intra tissue metabolite correlation plots, Pearson's linear correlation was calculated for data averaged at each time point for every pair of metabolites.
Using a stringent threshold (P<0.001), every metabolite in a given family that correlated with at least one metabolite in a target family contributed to ribbon width between the two families. Segment size reflects the number of detected metabolites in each tissue. Ribbons with width <10% of family size at both ends were omitted for clarity. Data was processed in MATLAB R2019a (Mathworks, US) and plots made using the Circlize visualisation package in R (12).

Statistics
Statistical tests and sample numbers are specified in figure legends where appropriate. Full details of statistical tests are provided in Supplementary Dataset S1. Statistical tests were conducted in GraphPad Prism. Throughout, * denotes p<0.05, ** denotes p<0.01, *** denotes p<0.001 and **** denotes p<0.0001. Plots were produced in GraphPad Prism, using the R package ggplot2 (13), or using Matlab.

SUPPLEMENTAL DISCUSSION AND STUDY LIMITATIONS
Characterising rhythmic processes is complex with numerous methods available, each with strengths and caveats. Here we use a comparative analytical approach (CompareRhythms, (3)), which has improved implementation of cosinor regression for comparing rhythmicity directly between conditions. Utilising cosinor analysis and model selection, as done within CompareRhythms, is a much needed development in differential rhythmicity analysis. JTK-cycle (14) and similar methods such as RAIN (15) are still commonly used for this kind of comparison.
However, these methods are not built for identifying differences in rhythmicity between conditions. JTK-cycle, RAIN and similar methods are built for rhythm detection based on samples from one condition. Using them for differential rhythm detection by comparing 'rhythmic' gene lists from each condition (analysed separately) leads to a high rate of false discovery (16); for example, just achieving significance in one sample, and just failing to reach significance in another, results in the incorrect inference that they are different from each other.
Model selection methods such as CompareRhythms (3) and dryR (17) classify paired gene profiles into distinct categories: arrhythmic, loss of rhythm, gain of rhythm, same rhythm or changed rhythm. The gene profiles are then modelled jointly across the two conditions and probability scores for membership of each category obtained. Whilst we have predominantly used this comparative analysis approach, we have also included analyses with more traditional methods (JTK cycle) to make our study fully comparable with previous work. Both methods estimate similar numbers of genes to be rhythmic overall, and identify similar changes in rhythmicity, functional pathway enrichment and potential upstream regulators. This is robust across different expression and probability thresholds.
Nevertheless, we must acknowledge that using different thresholds and/or analysis approaches will change the number of rhythmic genes detected within a tissue and given condition. Indeed, drawing conclusions based on absolute gene numbers assigned to any given category (rhythmic, gain, loss etc) should be avoided. Importantly, our gene ontology analyses and upstream regulator analyses were consistent across multiple methods, expression level cut-offs and probability thresholds.
In this study, we analyse rhythmic changes in a complex disease model which can show considerable differences between individual mice. Therefore, variability in disease state between mice could influence the assessment of differential expression and rhythmicity between naïve and CIA mice. To mitigate this, we implemented strict criteria regarding disease severity for sample inclusion ( Figure   1A), and include robust numbers of replicate samples collected over multiple independent experimental runs for each time point. Our transcriptomic and metabolomic samples were collected from four independent experimental mouse cohorts, and standardised to ensure disease severity was evenly distributed across time points and replicates. We characterised five samples at each time point to minimise the risk of false positive detection due to noise and biological variability (18).
We cannot rule out that behavioural effects associated with CIA (such as reduced locomotor activity) may have contributed to transcriptional and/or metabolic differences between naïve and CIA mice. Our telemetric assessment of activity and body temperature suggests that overall levels of activity are reduced in CIA mice once they develop symptomatic disease; however, there remains a significant difference in activity between the dark and light phases. Due to the method of activity measure (radio telemetry using DSI TA-F10 remotes), we cannot determine absolute activity levels in the animals (as the activity counts generated are not directly proportional to distance travelled). Importantly, we show that physiological measures (body temperature) and peripheral entraining signals (corticosterone) remain robustly rhythmic in these animals even during symptomatic disease.
Previous studies have found that exercise can entrain the mouse circadian clock (15), and can alter rhythmic physiology and clock gene expression in peripheral tissues (16), including skeletal muscle (17). In light of these findings, it is notable that we do not observe changes in rhythmicity of the core clock genes in muscle tissue with disease (Supplementary Figure S3C). We therefore consider it to be much more likely that the changes we observe in clock gene expression in the joint, the primary site of inflammation, are attributable to disease processes rather than any loss of activity-related entrainment.     *** * ** ** *** **** **** *** **** **** **** **** **** *** **** **** **** ****     Venn diagram demonstrating the overlap between genes differentially expressed in CIA mice and targets of STAT3 binding, as determined by chromatin immunoprecipitation (19). Cistrome data was extracted from CistromeDB (20).
STAT3 targets were defined as genes for which the cistromeDB score was greater than zero in at least one IL6-treated hepatocyte sample.